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Abstract 

Geometric integration theory can be employed when numerically solving ODEs or PDEs with con- 
straints. In this paper, we present several one-step algorithms of various orders for ODEs on a collection 
of spheres. To demonstrate the versatility of these algorithms, we present representative calculations for 
reduced free rigid body motion (a conservative ODE) and a discretization of micromagnetics (a dissi- 
pative PDE). We emphasize the role of isotropy in geometric integration and link numerical integration 
schemes to modern differential geometry through the use of partial connection forms; this theoretical 
framework generalizes moving frames and connections on principal bundles to manifolds with nonfree 
actions. 

1 Introduction 

In this article, we describe a set of algorithms on multiple copies of (with possibly nonlinear interactions), 
with applications in material science. The physical processes we are interested in are modeled by PDEs of 
the form 

d 

— /i(x,i) = A(/x(x,t)) X /2(x,t), /2(x,0) = /Xo(x). (1) 

Here /x : 23 ^ M'^ is the field describing the physical process of interest, with x e 23 denoting the spatial 
variable on some closed, compact subset of M.^, and A(/x) is a (typically nonlinear) function of /x. One im- 
mediately notices that ||/x(x, •)|| is constant for all time; this constraint should be respected by the numerical 
methods used to study such systems. Classical numerical integrators fail to preserve this norm, and this 
motivates our use of geometric integration. 

The paper is organized as follows: we begin by introducing a problem in micromagnetics. We then describe 
the Lie group framework suitable for the general problem ([T]). The new, and fairly geometric, mathematical 
constructs of partial moving frames and partial connections, which can be used to select generators for use 
in numerical methods, are introduced in section 3. Readers who are unfamiliar with modern differential 
geometry, particularly principal bundles, are advised to skip this section at a first reading. In section |4j 
we present some one-step algorithms of different orders. An arbitrary function appears in these algorithms; 
different choices of this function yield distinct discrete trajectories. We describe a possible choice that 
is related both to moving frames ([D [21 and, in the case of a geometric version of the forward Euler 
method, to discretization error minimization. Work is in progress to identify analogous function choices 
for higher order methods ( 4J). We present the results of numerical experiments carried out using some 
of these functions for the symplectically reduced free rigid body, a Hamiltonian system on the two-sphere. 
The conservative nature of the rigid body system facilitates both the derivation of appropriate functions 
and the assessment of the relative performance of the algorithms. Finally, we present some representative 
numerical results for micromagnetics that indicate that these geometric integration schemes are competitive 
with conventional numerical algorithms. Our numerical experiments were carried out using first, second, 
and fourth order methods. 

** Department of Mathematics, University of California, Santa Cruz, ^Department of Mathematics and Statistics, McGill 
University, Montreal 
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1.1 Why use geometric integrators? 

One of the intrinsic features of the system ([T]) is that the vector field /x(x, •) evolves on a sphere. However, 
a classical integrator will update /i(x, t) at time t using the approximation 

/x(x, t + At) w /x(x, t) + F(/i(x, t),t, At). 

The particular form of F(-, •, •) depends on the algorithm chosen; however, it is clear that such updates 
correspond to translations of /i(x, t), not rotations. Thus, a classical integrator does not account for the fact 
that fj, evolves on a sphere. This constraint is difficult to efficiently impose in practice. 

A naive approach is to keep track of changes in the norm ||/^(x, t)\\ during a numerical experiment, and 
renormalize the iterates after a prescribed tolerance has been exceeded. However, this renormalization is 
equivalent to the aphysical addition (or subtraction) of energy to the system and is therefore an undesirable 
solution. In addition, this renormalization would also affect A(/i). We shall see that in the context of 
micromagnetics, this change is nontrivial and nonlinear. 

Another observation is that the component of A(/i) parallel to ^ does not influence the solution curve 
/^t(x, t). However, this component does alter the discrete trajectories generated by numerical algorithms. An 
appropriate selection of the normal component can improve the performance of the scheme. Preliminary 
numerical investigations suggest that such corrections may allow increased accuracy in the capture of key 
features, e.g. orbits and/or conserved quantities, at a moderate computational cost, without the introduction 
of significant numerical artifacts, e.g. numerical accelerations. 

1.2 Who needs geometric integrators? 

The philosophy that numerical algorithms should respect properties intrinsic to the system — momentum 
conservation, evolution on a manifold, Hamiltonian structure — is not a novel idea. (See, e.g., [S], [5], [7]). 
Many classes of algorithms have been developed with this guiding principle, especially for energy- conserving 
or symplectic systems. As far as we know, the application of such methods in the context of numerical 
micromagnetics is relatively new. Encouraged by the success of these algorithms in an industrial context 
([5]), we believe that these techniques will find wider use in the material science community. 

Modern magnetic materials are used in an increasingly large number of applications, including thin film 
read heads and recording media nanocrystalline permanent magnets, and magnetohydrodynamic fluids. 
In addition, there has been much interest in the use of "smart materials" , including magnetostrictive ac- 
tuators and organic ferromagnets [10]. These magnetic materials exhibit different responses corresponding 
to varying magnetic fields. For example, the resistance of a read-head ferromagnetic sensor changes as the 
device rotates over a recording medium. The magnetization of the material directly interacts with the other 
physical and chemical characteristics of the material; micromagnetics theory describes this interaction at a 
microscopic level. Of particularly great interest is the correlation of physical microstructure and magnetiza- 
tion; the ability to predict the response of one to variations in the other is crucial to the further development 
of these materials and devices. 

Another physical system whose mathematical model strongly resembles that of micromagnetics arises in 
the study of liquid crystals. Nematic and smectic liquid crystals form the basis of the operation of many 
every-day devices, such as LCD's, telecommunication devices, thermometers, projection systems — even 
mood rings. The devices operate on the principle that a suitable applied field will change the orientation of 
the liquid crystals, while conserving the pointwise-norm. Therefore, the resulting mathematical model has 
the same constraints as those in the micromagnetics situation. This system has been studied extensively, 
yet a norm-conserving algorithm has only been presented for an extremely simple situation ([TP). 

The study of long-chain molecules such as those occurring in bio-molecular systems yield another possible 
application area for the methods we develop. These chemical systems have more complicated constraints 
on the geometry of possible configurations; the configuration sought is one that minimizes a free-energy 
functional. Statistical thermodynamics considerations permit the reformulation of these optimization prob- 
lems as evolutions in time on given manifolds to a steady state ([T^,[I3])- Such studies are used in the 
development of new drugs. 
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1.3 A mathematical model of micromagnetics 

At any point in time, the magnetization fi is constant in small regions in the material (termed a domain), 
and the switching of these domains from one state to another is the basis of the functioning of devices built 
with these magnetic substances jl4| . Ideally the entire device would be one large domain, which would 
switch instantaneously when an applied field was imposed. In practice, there are several domains, and 
the net magnetization in a desired direction is not optimal. We are interested in tracking the evolution of 
these domains, which entails following the local behavior of fi. In industrial applications, a ferromagnetic 
device is subjected to changing magnetic fields (corresponding to the various applications — a disk moving 
underneath a read- head, changes in applied voltages for magnetorheological fluids, etc). One is interested in 
the response of the device to the gradual changing of these external fields. 

A model that is widely used in the industry is the Landau-Lifshitz-Gilbert (LLG) model of micromagnetics, 
which describes the evolution of the state of magnetization /x in a ferromagnetic sensor, occupying a region 
25 in space. The LLG equation for the magnetization /x(a;, t) is given by 

d 

— /^t = -/X X Hcff(/x) - A/x X (/X X Hcff(/x)), ||/x(x)|| = 1 V X e S. (2) 

Here A is a damping parameter and H^ff is an effective magnetic field, described in detail below. The first 
term on the right hand side of ^ describes the (undamped) Larmor precession of /x about Hoff and is derived 
from first principles [151 116) . It is observed in physical experiments, however, that changes in magnetization 
decay in finite time. The second term in ([2]) is a phenomenological term (called the Gilbert damping term, 
see [ddH]), added to describe this damping behavior; it cannot be derived from first principles. There are 
situations under which this system is stiff (see, for example, [121 [20] ) , and issues of numerical stability of the 
integration scheme are therefore important. We also thank one of the referees for pointing out |21j . where it 
is suggested that poor representations of the exchange energy may lead to the observed stiff behavior. 

The effective field, which causes the magnetization to change, is derived from energy considerations [15] 
and varies nonlinearly with /x. More precisely, 

Heff(/ii) = AA/x + /Xq (-V0 + Happ) + K{fi ■ e)e. (3) 

The parameters A and K are material constants of the permalloy being studied, and /Xq is the permeability 
of free-space. The field AAfi is called the exchange field, preventing rapid spatial variations of /x and the 
formation of arbitrarily fine domains (see [221 123j for examples on how this term is computed in general) . The 
final contribution in (|3| is due to the nature of ferromagnetic crystals, which causes the magnetic moments 
to align in preferred directions. This effect is incorporated in the LLG model through the uniaxial anisotropy 
field K{fi ■ e)e. The external applied field is denoted by Happ. The nonlinear, nonlocal contributions of /x 
arise through the demagnetizing field, — V(/), where </> solves the Poisson problem with suitable boundary 

= V • /X in M^; [0] = 0, 

and radiation conditions 

(f> — o ^ I — j-^ as |x| oo. 

Here [u] denotes the jump of the function u across 923. Many different methods exist for the calculation of 
this field, including the use of the full Maxwell system, FFT techniques, finite element methods, multigrid 
approaches, finite differences, and recently, fast multipole methods, [2Il[^ [^[?ni ^ l ^ l5 ^ l5T 1 [5 ^ [55 1 |Ml 
[55] , among others. While there are still several unsettled issues in the area of demagnetizing field calculations 
(see [30j for a sharp critique on existing methods), it is not our intention in this project to duplicate this 
work. Instead, we shall focus on developing a time-stepping method that is robust, accurate, and requires 
relatively few field evaluations. 

Theoretical developments in micromagnetics are driven by industrial demands, and the need for accurate 
algorithms is now imperative. Conventional algorithms are still being employed for highly sensitive calcu- 
lations on large sensors, and are becoming increasingly inadequate. Moreover, the time-scales inherent in 



dn 



= /X • n on d'B, 



(4) 
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these problems vary from nano-seconds (in disk drives) to tens of seconds; hence integration techniques that 
remain effective over long times are required. An ideal integrator would resolve solutions accurately over 
very small time steps, while allowing large time steps to be taken when the system evolves over a long period. 



1.4 Potential problems with renormalising 

We have already described some of the issues with classical integrators with respect to norm conservation. 
Here we show that renormalizing /x to conserve the norm while using conventional integrators changes the 
potential (/) in a nonlinear fashion; it is easy to construct a simple example in which the renormalization 
introduces a significant change in the demagnetizing field. Assume that the magnetization /x satisfies 

Kx,y,z) = ai + b{x)j, (5) 

for some constant a and scalar function 6; here i,j denote the usual unit vectors in the x and y directions. 
Before renormalization, the field M is divergence-free: V • /x = 0. However, the rcnormalized field is not 
divergence-free: 

H \ abb' 



(a2 + 62)3/2 ■ 

The potentials obtained by solving Q are clearly not the same for the original and renormalized fields. The 
effect of renormalization on the demagnetization field is particularly significant near domain boundaries. 
For example, if the function b appearing in ([S]) is a step function, then the divergence of the renormalized 
magnetic field is a delta function. We thank F. Reitich, [36 , for this illustration of the dangers of normalizing 
vector fields. 



2 Lie group methods for the system jj, = A(/x) x ^ 

We recall that the general system under consideration is 
d 

— At(x,t) = A(/x(x,t)) X /x(x,t), /x(x,0) = /Xo(x), x e S. (6) 

The initial condition /Xq satisfies ||/Xq(x)|| = 1 for all x e 23, i.e. /Xq : 23 — > 52, the unit sphere in M"^ Since 
^||/x(x, i)|| = for all x and t, ^{ ,t) : 23 — > 52 for all t, the rotation group SO{3) acts transitively on 
S'2, i.e. any point on the sphere can be rotated onto any other point on the sphere; hence there are time 
dependent curves Q : 23 — 5*0(3) satisfying 

M(x,t) = Q(x,<)/Xo(x), VxeS. (7) 

We emphasize that these curves are not uniquely determined by ([7|. 

In the spatially discretized version of the system, we choose N grid points x„,n = f, 2, in 23 and 
consider only the values of the magnetic field at those grid points. Given a curve /Lt(-, t) : 23 — > 5*2, we define 

M{t) := (M(x„i))tieM:= (S^ . 

Where there is no confusion, we suppress the argument t. The fully discretized version of the system ^ is 

M = A(M) X M. (8) 

Here and throughout the paper, vector operations such as cross or inner products on (5*2)^ or (M^)^ should 
be understood as the usual operations in M"^ performed on each of the A^ component vectors. The Lie group 
S = (iS'0(3))^ acts transitively, but not freely, on M. That is, any point in M can be mapped onto any other 
point in M by the action (by component-wise rotations) of S on M, but the group element accomplishing 
any such transformation is not unique. The isotropy subgroup of a point M in M ( where M is the group of 
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transformations fixing M) is an A^-dimensional torus. As in ([t]), there are smooth time-dependent curves Q 
in the group S satisfying 

M(t) = Q(i)M(0), and hence M(t) = Q(t)M(0). (9) 

We can use the Lie algebra 

© = (so(3))^ = {skew symmetric 3x3 matrices}^ w {M?)^ 

of S, which is the tangent space to S at the identity, to put ([£]) into a more famihar and convenient form. 
The identification of so(3) with is implemented using the map skew : — > so(3) given by 

skew [^] := 

i.e. skew [^] x = £ x x for all x G M'^. The matrix commutator bracket on so(3) corresponds to the cross 
product on under this identification. 

Given any differentiable curve Q{t) in S, there exists a curve ^(t) e (M'^)^ satisfying 

Q(t) =skew [i{t)]Q{t), 

where the product of skew [^] and Q is the usual matrix product. Thus the system ([8| is equivalent to 

M(i) = Q(t)Mo = skew [i{t)] Q(t)Mo = £,{t) x M{t). (10) 




Comparing (10) to ([8|), we see that 

|(t) = a;(M(i)), where a;(M) = A(M) + cr(M)M (11) 

for an arbitrary scalar function cr : M — > M. The flexibility in the choice of map a arises from the non-freeness 
of the action of 5*0(3) on 5^, and thus the action of S on M; distinct ODEs 

Q — skew [a;((5 m)] Q and Q = skew [u3{Q m)] Q, 

where 01(771) — cD(m) e span[m] for all ra ^ S'^ , will typically have distinct solution curves in 50(3), but the 
images in 5^ of those solution curves under the map Q ^ Q ■ mo will coincide. 

When numerically simulating ([s]), we want a time-stepping method that ensures that M„ G M, i.e. that 
the norms of the component vectors are identically equal to one. We can regard (10) and (11) as defining a 
family of ODEs 

Q = skew [a;(gMo)]Q, 

parametrized by Mq G M and cr : M — t- M, on the group S and use the techniques developed for geometric 
integration on Lie groups to determine approximate discrete solution curves of these ODE (see [STl 1551 15^ HO] 
and the references therein). When combined with the action of 9 on M, these techniques yield geometric 
integration schemes for ^ that exactly preserve the constraint M„ G M, regardless of the step size or the 
order of the integrator. 

The key idea is the following. Suppose we are given a (right) trivialized form g — (,{g)g of an ODE on a 
Lie group G and an algorithmic exponential Exp : g — >■ G mapping the Lie algebra g of G into G. Then an 
integrator of order k corresponds to an update of the form 

gn+i = Exp(i^(5„_p, ...,gn, ^t))gn, 

for some map F : G^+^ x M — > g determined by the algorithm and the generator ^. Here we consider only one 
step methods, with p = 0. We emphasize that the algorithmic exponential need not be the true exponential 
of the Lie group, or even a good approximation to the true exponential; all that we require is that it map 
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algebra elements exactly into the group G and that the algorithmic exponential of the zero vector in g equal 
the identity element of G. For example, the Cayley transform, given by 

cay(l) = (/ + skew [4/2] )(/ - skew K/2])-\ 

is an algorithmic exponential for the rotation group that has long been used in computational mechanics. The 
Cayley transform has long been used to implement exact rotations in elasticity and plasticity simulations; 
see, for example, ^4TJ|35]). More recently, it has been utilized in the geometric integration of a wide variety of 
mechanical systems, including the LLG equations; see |37| [33 HSl EU US] , and references therein. Algorithms 
of arbitrarily high order can be constructed using the Cayley transform, despite the fact that is only a second 
order approximation of the matrix exponential of 50(3). (Note that the Cayley transform is, in fact, an 
algorithmic exponential for any matrix group determined by a quadratic constraint. See, e.g. [46].) For the 
rotation group 50(3) on , both the true matrix exponential and the Cayley transform can be efficiently 
evaluated and have frequently been used as algorithmic exponentials. The true exponential takes the form 

fe\ T , ^"^ll^ll 1 rti , l-cos||4|| 2 
exp 4) = / + -TTZT— skew I + -— skew ^ . 

Wiw Wiw 

The image of a vector x e M"^ under the action of cay(4) takes the simple, readily evaluated form 

^xx+i^x^xx) . 

Hence we shall use the Cayley transform (actually, N copies of the Cayley transform) as our algorithmic 
exponential Exp : © — )■ S- We observe that the Cayley transform has the advantage that the entries of cay(4) 
are rational functions of the components of 4; in particular, no trigonometric functions need be evaluated. 
To summarize this section, we have rewritten the discrete system 

M{t) = A(M) X M(t), M(0) = Mo 

as an ODE 

Q(t) =skew [u;(Q(i)Mo)]Q(t). 

on the Lie group S = (50(3))^. We now need to describe choices of infinitesimal update maps F and 
generators A that determine one step numerical updates Qn+i = Exp(F((5„Mo, At))(5„ and associated 
updates 

M„+i = Q„+iMo = Exp(F(M„, Ai))M„ 

with specified properties, e.g. a specified order of overall accuracy. The construction of suitable updates is 
the subject of the next section. 

3 Generator selection — Expansions, curvature, and partial con- 
nections 

A natural and obvious goal in selection of a numerical scheme is the achievement of the highest possible 
accuracy working within the given constraints. However, the prioritization of the constraints (efficiency, 
stability, developer effort, preservation of key features of the modeled system, etc.) can lead to significantly 
different approaches to the achievement of this goal and correspondingly different schemes. For the purposes 
of this discussion, we shall assume that we are given a family of one step methods of the form 

(M„+i), =cay(F(M„,Ai),)(M„), j = l,...,iV (12) 

for some map F : M x M ^ {M.^)^ . (Recall that M := (5^)^.) 



cay(4)x = x 



1 

1^ 
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We can optimize accuracy within this family of methods by selecting an appropriate isotropy algebra 
correction; specifically, we shall first compute the update generator given some 'default' choice of generator, 
and then use that generator to determine an element of the isotropy algebra of the current state that 
minimizes the discretization error for that update. We will derive conditions specifying this choice of isotropy 
element for algorithms on M utilizing the action of S using a traditional series expansion approach to the 
computation of the discretization error. (See [3] for the application of this approach and those discussed 
below to general homogeneous manifolds.) Subsequently, we will discuss some more geometric and, in some 
cases, less computationally intensive, approaches to this task. 

The geometric approach most closely related to the naive series expansion treatment is the use of geodesic 
curvature to characterize the essential information about curves on manifolds, e.g.; solution curves of differ- 
ential equations. A less directly related approach, but one that coincides with the direct error minimization 
approach for the forward Euler method on M is the use of a partial moving frame and the associated partial 



connection form to determine a choice of generator. We shall define these constructions in section 3.2 



3.1 Algorithms on 5^ 

To demonstrate the influence of the isotropy algebra on discrete trajectories, we first study the action of the 
rotation group 50(3) on a single sphere 5'^. The techniques we use in analysing the case of a single sphere 
immediately generalize to M An autonomous vector field X on S'^ satisfies {X{m),m) — for all m £ S'^; 
hence there exists a (nonunique) map A : S''^ — > M^, called a generator of X, satisfying 

X{m) ~ A{m) X m 

on S*^. Recall from our earlier discussion that distinct choices of generator A typically yield distinct discrete 



trajectories when used in a numerical algorithm of the form ( 12 ). We will compare the performance of various 
schemes with two different choices of generator. The first choice is the 'default' or 'natural' one, which is 
not assumed to have any particular geometric properties. For the rigid body system, we will take the body 
angular velocity as our default generator. For arbitrary vector fields X, there need not be a natural choice 
of generator; our intent here is to make a plausible choice of generator that one might make if the issue of 
isotropy were not taken into account. The second choice is the orthogonal generator, i.e. the unique map 
Ao-.S"^ ^ satisfying 

X{m) = Aoijn) X m and (ylo(TO),m) = 

for all m e S^. 

A general and direct, but potentially computationally intensive, approach to the choice of a generator 
that reduces the discretization error is to compute the lowest order nonzero term in the series expansion 
for the discretization error for the family of algorithms under consideration, leaving the isotropy algebra 
component of the generator as an undetermined parameter, and thus determine conditions on the isotropy 



component that minimize the error. See [3 for a general treatment of this approach and section 4.1 for the 
derivation of the optimal generator choice for our forward Euler algorithm for the LLG system. Here we 
briefly explore some alternatives to this approach that have natural geometric interpretations. 
We consider an Euler update of the form 

^At{m) = Exp(At (Aoim) + a{m) m)), (13) 

where cr determines the isotropy (normal component) contribution. In the special case that the algorithmic 



exponential is given by a rescaling of the usual matrix exponential, e.g. by the Cayley transform, ( 13 ) 
satisfies _ 

^^AtiiTT-) = exp(r(m, At) {Ao{m) + a{m) to)) 

for some rescaling r of time. Hence, in this case, 9^t(TO) is given by a rigid rotation of to about an axis 
depending only on to. Hence the curve 



,(to) ^ {jt(TO) : \t\ < e} 
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is a segment of a circle in S^. Our goal is to choose a so as to obtain the best circular approximation at m 
to the true orbit segment 

0,{m) ^ {9tim) : \t\ < e} . 

If X{'m) ^ 0, then the optimal circular approximation to Oc(m) at m can be characterized using the geodesic 
curvature 

{{X ■V)X{m),mx X{m)) 
~ ||X(m)||3 

of Oeim) at TO. The best circular approximation to ©^(to) at to is tangent to X(to) at to and has geodesic 
curvature equal to that of Oe(TO) at m. The first condition is clearly satisfied for any consistent update. The 
geodesic curvature fcg(m) of T^{m) is easily seen to satisfy |fcg(TO)| = | cot (f)\, where (f> is the angle between to 
and Ao{m) +a{Tn)m (see, e.g. [17], p. 249); thus |A:g(TO)| = |o-(to)|/||Ao(to)||. Hence optimal orbit capture 



within the class of updates ( 13 ) is obtained using 

tJcorM :=fcg(TO)||X(TO)||. (14) 

If Ot(m) is itself a segment of a circle, then (Tcor yields r^{m) — 0^{m). Hence any torsion-free orbits, 
e.g. the separatricies of the reduced rigid body equations, are captured exactly by this version of the Euler 
method. Note that the choice Aq is suboptimal for the Euler method unless kg = along the tractory of 
interest, i.e. unless the desired trajectory is a great circle. ^ 

For higher order methods, the axis of rotation used in the update map is typically time dependent and 
hence the corresponding algorithmic trajectory segment typically is not circular (i.e. it has nonzero torsion). 
Hence the simple argument used in the preceding paragraph cannot be applied. However, the strategy of 
curvature-matching can still be followed. Since a smooth curve on a two dimensional manifold in M"^ is 
determined up to a time reparametrization by its geodesic curvature, we can determine the conditions on 
the generator imposed by the restriction that the geodesic curvature of (™) match that of {'m) to some 
order. The higher order derivatives of the curvature can either be determined analytically for a given vector 
field X or numerically approximated using standard difference schemes. 

More generally, the choice a generator of a vector field on a homogeneous manifold can be viewed as 
a special case of the choice of a partial connection form, which generalizes to nonfree actions the classical 
connection form on a principal bundle. A partial connection form is a Lie algebra-valued one-form with 
appropriate equivariance properties. In section |3.2| we state the relevant definitions and present a family of 
partial connection forms on open subsets of S'^ that yield a discretization error-minimizing algorithm and 
determine an algorithm that captures orbits to second order for any dynamical system on a single copy of 
S^. The interested reader is referred to [3| for a more detailed treatment of partial connection forms and 
related constructions. 



3.2 Partial connection forms 

We now briefly discuss a general geometric approach to the selection of generators, using a generalization of 
the connection form on a principal bundle. For a more detailed treatment of this generalization and proofs 
of the assertions given below, see Lewis et al. [2002]. Let T be a principle bundle, that is, a manifold T acted 
on by a Lie group G. Let the action of g G G in J* be denoted hy g ■ p — ^g{p) — $p(g), "^g & G,p £ T. Let 

Q -.— TeG denote the algebra of G and let g ■ p := ^T^^p ■ ^ : ^ G |. Recall that a connection on a principal 
bundle T is a distribution T satisfying 

TpV^Q-p^Tp and Tp$g • Tp = Tg.p, 

for all p € CP and g £ G. Specification of a connection T is equivalent to specification of an equivariant 
g-valued one-form a, called the connection form, satisfying 

a o Te4p ^ id, i.e. a{p){£,y{p)) = ^ for aU ^ e g. 
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for all p £ y. By equivariance we mean that a o T<^g — Adg o a for all g G G. The connection T and 
connection form a are related by the condition ker[a(p)] — Tp for all p E y. (See, e.g., [IQ for a detailed 
presentation of the properties of connections and connection forms.) 

The equivariance properties of connections and connection forms typically cannot be preserved in the 
context of nonfree actions, hence we relax these conditions, requiring only equivariance with respect to 
specified representatives of the isotropy equivalence classes. A map /3:GxM— J-Gisa slip map if 
P{g, m) ■ m = g ■ m for all g £ G and m e M. A (singular) distribution H assigning a complement to g • m 
in T,„M to each point m G M is a partial connection if there is a slip map /3 satisfying 

for all g E G and m E M. A partial connection form with slip map /3 is a g-valued one-form a on M 
satisfying 

a(m) (77m(™)) = f? mod g™, i.e. re$m(a o 1;$^ - id) = 

and 

*^i(ff,m)"('^) = Ad^(g, „)a(m) mod 

for all g € G and m E M. 

One natural source of partial connections is a generalization of a moving frame, in the modern sense 
introduced by Fels and Olver ([1^, [10]), i-e- a smooth equivariant map p : CP ^ G on a manifold CP with a 
free group action. Recall that a principal bundle CP is trivial if there exists a global section, i.e. a smooth 
map E : M — >■ CP from the base manifold M intoCP such that each group orbit G • p in CP intersects I](M) 
exactly once and the projection tt : CP — ^ M satisfies tt o E = id. This condition corresponds to the existence 
of a moving frame. Specifically, a global section E and associated moving frame p are related by the equality 

E(7r(p)) = p{py^ -p 

for all p G CP. A global section determines a flat connection, namely the connection that assigns to any point 
E(to) the subspace T^E • T^M. If we introduce the notation D'^ : TCP — > g to denote the right trivialization 
of the linearization of a map 7 : CP — >■ G, i.e. 

D'^iSp) := Tp(i?^(p)-i o 7) Sp 

for any p G CP and 6p E TpT, then Dp is the connection form of the connection determined by E. 

The rotation group 5*0(3) acts transitively on S'^ and freely and transitively on the unit tangent bundle 
U{S'^) = {u G TS'^ : \\u\\ ^ 1}. The map p : [7(5^) SO{3) taking u E U^S'^ to the orthogonal matrix 
with columns (m, u, m x u) is a (left) moving frame with associated connection form 

Dp{5u) = m X Sin + (m x 5u, m) m, 

where 5u E TuU{S'^), with m — 7t{u) and Sm ~ TuTt Su. (Here tt : U{S'^) — ?> S"^ denotes the canonical 
projection.) Note that we will regard u both as a tangent vector to the sphere at m and as a unit vector in 

Moving frames can be extended to manifolds with nonfree actions as follows: A (smooth) map cj) : M. G 
is a (left) partial moving frame if 

0s(™) := Hg ■ m){(f){m))~'^ (15) 

satisfies 

4>g{m) ■ m — g ■ m 

for all g E G and m G M. A partial moving frame on a submanifold S of a manifold M with a G action 



is a map (j) : § ^ G satisfying (15) for any m G S and any g E G such that g ■ m E §. The trivialized 
linearization D4> of a partial moving frame : M ^> G is a partial connection form, with associated slip map 
P{g,m) — (j)g{m). We refer to the trivialized linearization D(j) as the partial connection form associated to 
the partial moving frame 4>. 
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If a group G acts transitively on a manifold M, then every group orbit is equal to the entire manifold. 
In this situation, a (partial) connection form a assigns to each tangent vector a generator of that vector, so 
that 

{a{m){5m))-]^(jn) = 5m 

for all 5m € T^M and all m € M. In particular, given a vector field X on M, the map a; := Lxct : M g, 
i.e. uj{m) = a{m){X(m)), satisfies 

w(to)m(™) = X{m) 

for all m G M. Hence (partial) connection forms can be used to construct geometric integration schemes 
on manifolds with transitive actions. We are currently investigating the role of geometrically motivated 
choices of partial connection forms in the design of efhcient geometric integration algorithms. As can be seen 
in the example discussed below, simple, natural choices of partial connection forms can lead to significant 
improvement in numerical performance. 

To illustrate the somewhat abstract geometric constructions described above, we now present a moving 
frame associated to the action of the rotation group 5*0(3) on the unit tangent bundle U{S'^) of the sphere 
S'^ and an associated family of partial moving frames on S^. The partial connection form (which, in this 
case, is simply a map from the sphere to M.^) of one of these partial moving frames yields the discretization 



error-minimizing generators used in the versions of the forward Euler method described in sections 4.1 and 
[Sj As we shall see, this partial connection form and the associated generators can be derived without the 
use of the expansion of the discretization error. 

Any unit vector field y on a submanifold M of S'^ determines a partial moving frame cj) — p o Y on M, 
with partial connection form 

D(f>{5m) = m X 5m, + {Y(m) x [DY{m) ■ 5m), m) m. (16) 

The map (j)g associated to g G 5*0(3) is (f)g{m) — g exjp{9{g, m) m), where 9{g, m) denotes the angle between 
g~^Y{gm) and Y{m). Given an ODE m — X{m), we can set Y{m) = X{m)/ ||X(r7i)|| on some set M C 5^ 
containing no equilibria (i.e. zeroes of X); in this case, ^16^ takes the form 



A,/ N/c N c {(5m ■V)X,m X X(m)) 

D(j){m){5m) ^ m x 5m + — ^ ^—^ m. (17) 

II^MII 

In particular, 

D<P{m){X{m)) = mx X{m) + kg{m) \\X{m)\\ m, (18) 
where kg{m) denotes the geodesic curvature of the curve m(t) in 5^. 



Following (18), we set 



The partial connection form ( 16 ) can be used to select the isotropy correction map a used in (131 



.^m^(T^-i kgim)ix)\\X{m){x)\\ X{m){x)^0 

<j[m)[x) .- <i ^ X{m){x)=Q ■ ^^^> 



In our numerical implementation (19), we approximate kg{m) using the identity 

(to, to X to) (cj, uj X m) 

ll™ll IpII 

for a curve to {t) in 52 with nonzero velocity m = uj x m, where uj is orthogonal to m, and replacing to and 
TO with finite difference approximations. 

4 A collection of geometric integrators 

Recall that an update of the form 

Q„+i = Exp(F(g„Mo, At))Qn, 
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where Exp : © — > S satisfies (Exp(|))j ~ cay(^j), determines a one step method on the Lie group 9, where 
F : S X K — > © is determined by the generator and the selected scheme, and At denotes the time step. Given 
an ODE M(t) — A(M) x M on M, we wiU construct updates of the form 

M„+i = Exp(F(M„, A<)) • M„, (20) 

It is our objective to identify several classes of infinitesimal update maps F, leading to algorithms of first, 
second, and fourth order. 



4.1 Discretization error reduction through choice of a 

As was previously discussed, the normal component of the generator A does not influence the solution curves 
of the original ODE. Thus, if we have a numerical algorithm of order n, this component does not affect the 
solution up to order n. However, it typically does appear in the higher order terms of the approximation, 
and theoretically a suitable choice of this component will reduce the discretization error. For the forward 
Euler scheme an optimal choice of ct, in the sense that this choice minimizes the discretization error, also has 
a natural geometric interpretation. Here we derive this map a using a direct discretization error calculation; 
in the following section we shall discuss various geometric considerations that can be used in the selection 
of the generator to be used in a Lie group integration scheme. 

We consider consistent algorithms using standard methods on the tangential component of The 
normal component is treated as a function of the tangential one; we shall see that a component of the local 
discretization error at second order can be eliminated by a suitable choice of the normal component of the 
lowest order term in ^. For the sake of simplicity, we consider here only the lowest order case, in which the 
discretization error of a first order method is reduced by an appropriate selection of a. This is a particular 
example of a more general result covering a large class of manifolds and higher order algorithms. (See [2111].) 
Work is in progress (Lewis, Nigam, Olver) to possibly extend these or related results to an even larger class 
of systems, including the full discretized LLG system. 

We begin by examining the flow Jt of the ODE M ~ A(M) x M on M. This flow satisfies 

3'At(M) = M + AiAxM+ — (^Ax(AxM) + AxMj +0{At^). 

If the algorithmic update J'a* : M — > M is given by 

JAt(M) Exp(F(M, At)) • M 

for some map F(M, Ai) := V A^-^ -(M) and Exp : M"^ w so(3) 50(3) agrees with the exponential map 
to second order (e.g. Exp is the Cayley transform), then 

3F^t(M) = ^/ + Aiskew[A] + iAi2skew[A]V0(Ai3)^ M 

= + At skew [^i] + ^ At^ (^skcw [i^] + skew [i^f^ + 0{At^)^ M 

= M + At^i X M+ — (^1 X (^1 X M) +^2 X M) + 0(At3). 

We now derive conditions on the terms ^j^, 42v--iii the series expansion of F yielding algorithms of 
increasingly high order. The consistency condition for S'a* is ]Pm(€i — A) — 0, where Pm denotes component- 
wise projection onto the orthogonal complements of the component vectors of M, i.e. ((PmOjj Mj) = 0; 
i — 1^ . . . ,N . If 9^At is consistent, then, setting ai := M), the local discretization error is 



JAt(M)- JAt(M) _ At 

{-aiA + ($2 - A) X M)) + 0(At2) 



- 2 (-^iMx(^ixM) + (^2-A)xM))+0(At2) 
At 
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The algorithm is thus second-order accurate iff 



(|2-A,A)=0 and ai(A, A) = - A) x M, A). (21) 

In our geometric version of the forward Euler method with F(M, At) ~ A(M), = for j > 1; thus this 
method will not be second order. However, we are free to choose txi so as to satisfy the second equality in 



(211, e.g. 

(AA(M, Ai), A(M) X M) , , 
||A(M)||' 

where AA(M, Af) is some first order approximation to A(M) (e.g., a discrete difference approximation), 
yielding a discretization error-minimizing member of the family of algorithms with F(M. At) ~ A(M) + 
cr(M,Ai)M. 

Analogous expansions can be used to minimize the discretization error of higher order methods. However, 
the symbolic calculation of such expansions for high order schemes is, at present, relatively laborious and 
does not seem tractable for systems such as the LLG equations, in which the generator A is determined in 
part (the demagnetization field) by a nonlinear PDE. 

4.2 First order methods 



Using ( 20 1 , we now define geometric one-step methods that are natural analogs of the standard explicit and 



implicit Euler methods: 

f A(M) + cr(M, At)M forward Euler, 

F(M,Ai)=:<^ ' ^ ' (23) 

I A(M, At) + ct(M, At)m implicit Euler, 

where A(M, At) denotes the solution of the implicit equation ^ = A(Exp(At ^)M) and the scalar functions 
a and a are as yet unspecified. 

The numerical results presented in sections [3] and |6] illustrate the effect of the parameter a on the 
discrete trajectories determined by the forward Euler algorithm when applied to rigid body dynamics and 



the LLG micromagnetism model. We shall see that in the rigid body system, a satisfying ( 22 ) yields second 
order accuracy in energy tracking, and thus second order orbit capture for this conservative system. In the 
micromagnetics simulations, where damping plays a crucial role in the long term dynamics, large values of a 
cause the trajectories to sharply diverge from those of the ordinary forward Euler; however, the final state is 
the same. A closer look at the LLG equation shows that a larger value of a corresponds to the inclusion of 
more precession in the trajectory. These numerical results clearly show that different choices of the parameter 
cr lead to significantly different numerical trajectories and thus motivate the search for an "optimal" value 
of (T. In section [3] we describe a general geometric approach to selecting values for cr; in section ["4.1| we show 
that when used with the forward Euler method, this choice of a minimizes the discretization error. (See [3] 
for a description of this approach for more general manifolds.) 

4.3 Second order methods 

We consider four second order methods modeled on the classic Heun (RK2) algorithm. 
In the first, we use the 'default' generator A in the Heun method, i.e. 

Ff3f2(M,At) := i(A(cay(AtA(M))M) + A(M)). 

The second method, F^.^-^, is entirely analogous, but with A replaced with the orthogonal generator A^, 
where 

Ao(M) := A(M) - (M • A(M))M, 
and hence ((Ao(M))j , M^) = 0, j = 1, . . . , iV. 
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In the third method, the infinitesimal rotation determined by applying the Heun method to the default 
generator A is modified by addition of an appropriate isotropy element to yield a higher order of orbit 
capture; specifically, 

Ff,^2(M, At) Fj^2(M, At) + At' ^def (M)M 

The fourth algorithm is analogous, but with F^J^^ replaced by F^^.^^ and ddcf replaced by an appropriate 
function cTorth- (The function tTorth is a rational function in m and the components of A, but is significantly 
more complicated than (Xdof ■) Note that isotropy plays a role both in the choice of the generator and in the 
selection of a correction term. 



4.4 Fourth order methods 

In this subsection, we describe two families of fourth order Lie group integrators on M. We emphasize 
that these algorithms map a point M„ £ M exactly into M; they are fourth order accurate in the sense 
that they approximate the true trajectories within M to fourth order. A direct implementation of the 
classic fourth-order Runge-Kutta method to a vector field on M fails to maps exactly into the manifold M, 
while application of the classical RK4 method to the generator of the flow, followed by application of the 
exponential map and the group action necessarily yields an update in M, but typically does not give a fourth 
order approximation of the true flow. The generator must be modifled to account for the trivialization of 
the tangent bundle of the group; this modification can be implemented either before or after the stages of 
the Runge-Kutta method are computed and averaged. The first method we use is the RKMK4 method, a 
Runge-Kutta style method due to Munthe-Kaas ([39], [40]) in which each stage of a traditional RK4 method 
is modified so that the resulting generator, followed by (algorithmic) exponentiation and application of the 
group action to the manifold, yields a fourth order method. The second method utilizes a series expansion 
of the generator along the true flow, followed by a single modification to account for the trivialization of the 
tangent bundle of the group S, again followed by (algorithmic) exponentiation and application of the group 
action. 

We implemented the RKMK4 method [39], using the Cayley transform rather than the true matrix 
exponential. To implement a Lie group integrator for ^ using the Cayley transform, we make use of the 
fact that, for sufficiently small t, there is a function f : M ^ M'^ satisfying 

M(t) = cay(f(t))M(0). (24) 



Differentiating (24) with respect to t, we obtain 



M{t) = dcayf(f'(t)) X cay(f(t))M(0) 
= dcayf (f'(t)) X M{t) 
= A(M(<)) X M(t), 

where the map dcayj = I?cay(f) : — > M'^ is the right trivialization of the tangent map of the Cayley 
transform. 

Hence f ' and A are related by 

dcayf (f'(t)) X M{t) = A(M{t)) x M(i), 

which is equivalent to 

f'{t) = dcay^i(A(M(t)) + a{t)M{t)) (25) 



for some function a. The initial condition for (24) is f(0) = 0. 
The map dcayj : — >■ satisfies 

dcay, :=— ^(/+^skew[f]) 
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and 

dcay^^ ^I - - skew [f] + -f f"^. 



2 ' ' A 



The Cayley version of the RKMK4 method is essentially the conventional RK4 applied to ( 25 ) : having found 
M„ at time t„, we construct the update 

M„+i = M(t„ + M) = cay (F^^^^^^^^ ^^^-j 

where 

F«^4(M„^ At) = ^(Fi + 2F2 + 2F3 + F4) 


and 

Ai=AtA(M„), Fi =dcayo'(Ai) 

A2 = At A(cay (iAi) M„), F2 - dcayT^^CAa) 

A3 = At A(cay (^Aa) M„), F3 = dcayi^^(A3) 

A4 = At A(cay (A3) M„), F4 = dcayA3(A4) 

For more details on this method, see [55] . 

The second method involves a series expansion with respect to time of the generator A(5't(M)), where 
3^f denotes the exact flow at time t. Given the generator A, this expansion is computed by iteratively sym- 
bolically differentiating A(3'((M)) and then substituting A(M) x M for M. The third order approximation 
to A(5't(M)) is then modified to take into account the trivialization of the tangent bundle of 5*0(3) and the 
action of 50(3) on 5^. The specific expressions for this modification for the rigid body equations on the 
sphere are given in f|5j Given the implicit and highly nonlinear nature of the LLG equations, symbolic cal- 
culation of the derivatives of A for this system seemed excessively complicated; hence we did not implement 
this algorithm for the LLG system. 



5 An example: the rigid body flow on a sphere 

We now apply the results outlined above to a simple and familiar system, the reduced rigid body equations 
on the sphere. Given a positive definite symmetric three by three matrix I, define the vector field 

X{m) = 777, X 1^^777 (26) 

on 5^. This is a Hamiltonian system with respect to the Kostant-Kirillov-Souraiu symplectic structure 

51(777) (^ X 777,77 X — (™! C ^ 1) 

and Hamiltonian 

i7 (777) = ^(777,1-^777). (27) 



The system (26) is the symplectic reduction of the free rigid body equations on r*50(3); more concretely, it 
is the restriction of Euler's equation for the body angular momentum to the unit sphere. (Since the norm of 
the body momentum is preserved by the dynamics of Euler's equation, all spheres centered at the origin are 
invariant submanifolds.) The conservative nature of this system makes it particularly easy to measure the 
error in orbit capture; if the body is triaxial, i.e. the eigenvalues /i, I2, I3 of the inertia tensor I are distinct. 



the level sets of the Hamiltonian (27) exactly determine the orbits of the system. Thus in this situation 
the error in the orbit is a function of the fluctuation in the energy. As the numerical results given in tables 
1-4 demonstrate, geometric integration techniques yield efficient, accurate orbit capture for the reduced free 
rigid body, with good performance even for very large time steps. Note that the same randomly generated 
initial conditions and inertia tensors are used in all of the numerical simulations. 
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If the rigid body is axisymmetric, then all true trajectories consist either of equilibria (the 'poles' and 
the 'equator') or of steady rotations in the plane of symmetry. Note that in this situation, even an exactly 
energy-preserving scheme may allow drift across the family of one-point orbits along the equator. However, 
all of the methods considered here detect equilibria as such. Thus even in the axisymmetric case, we can use 
the energy to monitor orbit capture. We shall see that for some of the algorithms considered here, there are 
significant differences in performance on triaxial and axisymmetric bodies. Symmetries play a crucial role in 
algorithm design and analysis; see, e.g. [43j . However, we shall not explore those issues in any detail here. 



5.1 Euler methods for rigid body dynamics 

We now consider implementations of the families of algorithms described in f|4]for the rigid body equations. 
We take as our default generator the body angular velocity A{m) — l^^m. We first consider three first order 
methods, with infinitesimal updates 



m 



i^^"^(m) = I — {rn, I ^m) m — A{ni) — 2 H{m) m 



{X{m)J-^X{m)) 1 , r{u(m)) 

m = 11 TO H -—— — — TO, 

T(llu(TOjj 



where 



mmw 

and M : 5*^ — > M"^ are given with respect to an eigenbasis of I by 
r(x) = xi + X2 + and u{m)i := (Ij — Ik)'^Iimjmk 

for any cyclic permutation k) of (1, 2, 3). 
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Table 1: Maximum energy error over the trajectories given in figure [5?T 
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Table 2: Average global energy errors over ten sample runs with randomly generated initial conditions and 
inertia tensors, integrated over the interval [0, 100] using versions of the forward Euler method. 

In table [2] we provide the average maximum errors in the energy for time steps At — 10, 1, .1, and .01, 
using for ten randomly generated initial conditions and inertia tensors each for triaxial and axisymmetric 
bodies. 

The separatrix is exactly captured if the infinitesimal updates -F^th ^rar'; which coincide on the 
separatrix, are used. On the other hand, when was used to integrate ten sample trajectories with 
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Figure 1: Sample trajectories computed over the iutcaval [0, 200] using the time step At = 0.1 and, left 

to right, the first order infinitesimal updates Ff^^, -F^th' ^^"^ -^mr'- The upper row is computed using the 
inertia tensor of a trixial rigid body, while the lower row is computed for an axisymmetric rigid body. 
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initial conditions at random points on the separatrices of rigid bodies with randomly generated inertia 
tensors, the average errors over the integration interval [0, 500] were: 9.72 10~^ for At = 1, 3.89 10~^ for 
At = .1, and 1.94 lO'^ for At = .01. 

In the axisymmetric case, the forward Euler method with infinitesimal update associated to second 
order orbit approximation yields the exact solution when the true exponential map is used as the algorithmic 
exponential. (If the Cayley transform is used as the algorithmic exponential, then the orbits are captured 
exactly, but the algorithmic trajectories differ from the true trajectories by a time reparametrization.) Note 
that the 'default' generator and the orthogonal generator yield only first order orbit approximations even in 
the axisymmetric case. 

As implemented in our Mathematica code, the version of the forward Euler method with orthogonal 
algorithmic velocity is approximately 10% slower than the naive version, while the version that captures 
orbits to second order is approximately 30% slower than the naive version. 

5.2 Higher order methods for rigid body dynamics 
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Table 3: Average global energy errors over ten sample runs with randomly generated initial conditions and 
inertia tensors, integrated over the interval [0, 100] using versions of the Heun method. 



We implemented the four different versions of the Heun method given in section 4.3 for the rigid body 
system. Although the Heun methods F^^"^ and F^^^ described here are only second order accurate, the 
(local) discretization error in the energy is fourth order in the time step. 

The isotropy corrected versions used here take the form 



where 



^gl^{m, At) := Fg^im, At) + At^ adet(™)m. 



adef (m) := j^;^ , with <j J," := ~Ij {h + h){h - h 



for any cyclic permutation (j, k,£) of (1, 2, 3) and 

F„^f/(m, At) := ^.^^^(m. At) + At^ ao,.th(m)TO, 

where CTorth is another, significantly more complicated, rational function in m and the components of the 
inertia tensor. The isotropy corrections CTdof and Uorth given above determine algorithms yielding fourth order 
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Table 4: Average global energy errors over ten sample runs with randomly generated initial conditions and 
inertia tensors, integrated over the interval [0, 100] using several fourth order methods. 

energy capture. If the body is axisymmetric, -F^^^ preserves the energy to fifth order. As Table [s] shows, 
some of these algorithms appear to have better global energy capture than the single step discretization 
energy error analysis (which we carried out symbolically using Mathematica) would suggest. Plots of the 
energy errors in the sample integrations, with randomly generated initial conditions and inertia tensors, 
show that the energy oscillates about a very slow drift away from the correct value. Note that the energy 
correction term for the fourth order symbolic expansion method using the orthogonal generator is identically 
zero if the body is axisymmetric; hence the results generated by F^^^^^ and i^ocor coincide in this case. 

We consider six fourth order geometric methods. Four utilize a series expansion for the generator along a 
solution curve, while the other two use the RKMK4 algorithm (with the Cayley transform as the algorithmic 
exponential). Using the Cayley transform, the map F^^f determined by the default generator for the rigid 
body system on 5*2 is given by 




where A^i^m) = §jA{Jt{m))\t^o- The corresponding algorithm for the rigid body using the orthogonal 
generator is 

K^iirn,At) = ^^4^-i)(m) 

+ ^ (\\Aoim)fAoim) + (Jio(m), to) 
+ — (^Aoim),Ao{m)'^ Ao{m). 

The infinitesimal updates F^f and F^^^^^ can be modified by the addition of an appropriate multiple of the 
argument m to yield an additional order of energy, and hence orbit, capture. The scalar correction functions, 
which are rational functions of the components of m and the inertia tensor, were determined by symbolic 
calculation. 
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6 Application of geometric integration to numerical micromag- 
netics 



In this section, the geometric integrators developed in section |4] are applied to the Landau-Lifshitz-Gilbert 
equations of micromagnetics. The exact solution of this system is typically not available to us; indeed, it 
is the lack of precise analytical results for comparison which makes numerical micromagnetics a challenging 
field. In our examples, we chose the largest possible time steps for a given method that would lead the 
system to the solution computed by a higher order method (within prescribed tolerance). 

As mentioned earlier, numerical micromagnetics has attracted much attention in the mathematical com- 
munity, for several reasons. In this paper we are focussing on the time-stepping aspect of the problem. The 
application of geometric integration techniques in this context is relatively new, see for example, |44[ I45| . 
Recently, another technique which modifies existing integrators was developed for numerical micromagnetics, 
[5T] . This new integrator is of the "step-and-project" class, but is stable. 



Figure 2: The trajectories followed by the usual forward Euler At = 0.0001 and with the geometric forward Euler 
(/S.t = 0.01^ with optimal a are almost identical; if we assume cr = 10 (an arbitrary choice), the trajectory precesses 
more before reaching the final point. The figure on the right shows the evolution of the optimal a. The applied field 
is uniform and weak, specifically, Happ = (0.05,0.05,0). 

We were particularly interested in the behavior of the free parameter cr which appears in the geometric 



time-stepping algorithm (11). In the rigid body case the parameter can be chosen to improve energy con- 
servation. Here the system is dissipative and a criterion for the selection of a is not immediately obvious. 
For the forward Euler implementation, we can derive a relatively simple expression for a function a that 
minimizes the discretization error. For higher order methods, analogous functions can be described in terms 
of the series expansions of the true and algorithmic flows, but the cost of computing these expansions, partic- 
ularly for systems such as the LLG equations, rapidly becomes prohibitive. Work is in progress to determine 
computationally tractable criteria for the selection of the isotropy component for higher order methods. 

6.1 Description of the model problem 

We describe a model problem for the LLG, for which the analytical solution was particularly simple. Recall 
that the LLG for the magnetization fi(x,t) is given by ([2]), which we recall here for convenience: 

^M = -/^x Heff(M)- Vx (/XX Heff(/x)), ||/x(x)|| = 1 Vxe S. (28) 

Heff (m) = AAfi + ^0 (-V(/. + Happ) + Kifi ■ e)e. (29) 
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Figure 3: The evolution of one point in the ferromagnetic sample. The trajectories followed by the usual forward 
Euler method with time step At = 0.0001 and by the geometric forward Euler method with optimal a and time step 
At = 0.01 are almost identical; if we set cr = 10 (an arbitrary choice) in the geometric forward Euler method, the 
trajectory precesses more before reaching the final point. The figure on the right shows the evolution of the optimal a. 
The applied field is uniform, Happ ~ (5,0,0). 



In the experiments that follow, we set hq — K — A — 1.0 and vary the applied field. These parameter 
values are not taken from actual physical data, and were selected solely for purposes of illustration. The 
saturation magnetization was ||/x(x)|| = 1. 

We wish to construct a one-dimensional example where the computation of the demagnetizing energy 
would be simple. To this end, we assume the sample is contained in the infinite slab {{x,y,z) |0 < a; < 
l,y,z, G K^}. We assume the magnetization /x = M(a;), ie., the only variation in the magnetization is along 
the x-direction. Therefore, V • /x = (J^/^i,0, 0). We assume that there are 100 individual spins uniformly 
distributed along x G [0, 1]. These spins interact with each other through the exchange and demagnetizing 
fields. 

This example is admittedly a simplistic one; the true equilibrium solution for it can easily be found 
using analytical techniques. Therefore, the stopping criterion used was a comparison with the exact final 
equilibrium point. We see that the geometric integrators take trajectories which respect the point-wise 
constraints on the magnetization; conventional integrators do not. Thus, the paths traversed by these 
integrators will be different, as is seen in figures [2] and [3j As the step-size is shrunk more and more, the 
trajectories will converge. 



6.2 A first order method for the LLG 

In the first set of numerical experiments, we implemented the geometrical analog of the forward Euler algo- 
rithm for the LLG equation. We then tracked the evolution of the parameter a given by (22), an expression 
derived through arguments of discretization error minimization. We approximated the acceleration M„ of 
trajectories M(i) of ([8| using a one-sided discrete approximation of the derivative of M„ = A(M„) x M„. 
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Figures [2] and [3] describe the trajectories followed by the over-damped LLG system (without the Larmor 
precession term) for two different apphed fields. To address the issue of overall computational expense, we 
ran both a geometric forward Euler algorithm and the conventional forward Eulcr algorithm on the example 
introduced above. For each algorithm, we decreased the time step At until the trajectories converged to 
within 5%. We also required that the final equilibrium point corresponded to that computed by a fourth- 
order Runge-Kutta method with time step At = 0.0001, to within a relative error of 1%. The geometric 
forward Euler method yielded trajectories which converged, for this example, with time steps of At = 0.01 
and a CPU time of 1.77 seconds. The usual forward Euler required a time step At = 0.0001, with a CPU 
time of 3.88 seconds, to get similar behaviour. In addition, while the pointwise norm of M was conserved to 
machine accuracy by the geometric integrator, the standard forward Euler algorithm caused ||M|| to increase 
to 1.001183806 times its usual value by the end of the run. As a consequence, the trajectories traversed by 
the geometric and the usual algorithms differed, though they ended at the nearly same place. As we are 
only interested in the final equilibrium state of the system, we see the obvious merit of using the geometric 
integrator — we can obtain accurate final states while using much larger time steps. 

We see the effect of varying the scalar functions a on the trajectories is that of changing the amount of 
precession in the trajectory. We notice certain trends in the optimally chosen function a{t) in figures [2] and 
[3j and we shall investigate the relationship of these trends to the physical processes occurring at the same 
times in future work. 



Figure 4: An example of the full LLG system, with uniform applied field Happ = (5,0,0) and damping parameter 
X = 0.05. The trajectories were computed using the usual forward Euler method with time step At — 0.0001 and the 
geometric forward Euler method with ttme step At — 0.01. Here the trajectory taken by the geometric method differs 
appreciably from those of the usual method, though the final states appear to be similar. The lefthand plot shows the 
evolution of one point in the ferromagnetic sample; the righthand plot shows the evolution of the optimal sigma at 
that point. 

In figure [4j we implemented the code for the full LLG system, including the Larmor precession. The 
applied field is uniform, Happ ~ (5,0,0). The damping parameter A was set to a low value, specifically 
A = 0.05. The trajectories followed by the usual forward Euler and At = 0.0001 and by the geometric 
forward Euler with At = 0.01 and optimal a diverge appreciably, yet end at the same final state. The drift 
of the norm is now clearly visible (see figure ([5|). The usual forward Euler trajectory moves off the unit 
sphere in the standard Euler integration, while the geometrically integrated one does not. We see that the 
optimal cr now varies more (figure Hp). 
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Figure 5: Norm of magnetization in Fig(3). The geometric integrator exactly preserves the norm, even with a time 
step o/O.Ol. The usual forward Euler method shows a drift in norm, even with a time step o/ 0.0001. 

6.3 A second order method for micromagnetics 

In the next set of experiments, we implemented the geometrical versions of the Heim algorithm, derived 
in section |4j We did not have an analytical expression for the optimal choice of a. Therefore, we ran the 
experiments for several constant values of this parameter, and computed the order of convergence of the 
algorithm in At. 

The results were interesting, and rather striking. As a is varied, the order of convergence changes for 
the naive choice of generator. What should be noted is that the geometric algorithm appears to converge 
more rapidly than a conventional Heun method; the order of convergence was 0{At'^~^^), as was borne out 
in repeated experiments. 

The norm of M is conserved to machine precision for both generators. 

6.4 A fourth order method for micromagnetics: RKMK4 

We now present experiments with a fourth order method derived in section [4] Lacking an analytical expres- 
sion for the optimal choice of a, we ran the experiments for varying constant values of this parameter, and 
computed the order of convergence of the algorithm in At. 

In figurejsjwe track ||M|| over [0, 1] with a time step of 0.01. The classical RK4 method without projection 
shows a drift in the norm; this drift is of the order of 10~^, i.e., 0{At^). The Lie group integrator, on the 
other hand, shows no drift (up to machine precision). 

As we vary a, we observe that the rate of convergence of the algorithm varies, see (figure[9|. Again, there 
is clearly some optimal value of this parameter. This behavior is even more pronounced for the RKMK4 
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Figure 6: Convergence orders of the Heun method for varying a. We show experiments corresponding to two different 
generators. 

method than for the Heun method. 
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